- , vx ■ ^ J ' . 

■Mr'-.- T r.fai -’’.i- L ; .. 

ft. .In;'-; *V : \. .. - :;-. 

* ' ■ ’. " ;: -■ 



: i -m ^ 


<C 

on 


■■ - '- - 


CONTRACTOR 



NASA CR-1064 






GPO PRICE $ 
CFSTi PRICE(S) $ 

Hard copy (HC) . 
Microfiche (MF). 

ff 653 July 65 


AN ANALYSIS OF THE COUPLED 
CHEMICALLY REACTING BOUNDARY LAYER 
AND CHARRING ABLATOR 



f 

T f 

{ 


Part V 

A General Approach to the Thermochemical 
Solution of Mixed Equilibrium-Nonequilibrium, 
Homogeneous or Heterogeneous Systems 


by Robert M. Kendall 


Prepared by 

ITEK CORPORATION, VIDYA DIVISION 

Palo Alto, Calif. 

for Manned Spacecraft Center 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • JUNE 1968 


O 

LL. 

e 


(ACCESSION NUMBER 


(PAGES) 


(THRU) 

/ 


(CODE) 


(NASA CR OR TMX OR AD NUMBER) 


(CATEGORY) 



NASA CR-1064 


AN ANALYSIS OF THE COUPLED CHEMICALLY REACTING 
BOUNDARY LAYER AND CHARRING ABLATOR 

Part V 

A General Approach to the Thermochemical Solution 
of Mixed Equilibrium -Nonequilibrium 
Homogeneous or Heterogeneous Systems 

By Robert M. Kendall 

Distribution of this report is provided in the interest of 
information exchange. Responsibility for the contents 
resides in the author or organization that prepared it. 

Issued by Originator as Aerotherm Report No. 66-7, Part V 

Prepared under Contract No. NAS 9-4599 by 
ITEK CORPORATION, VIDYA DIVISION 
Palo Alto, Calif. 

for Manned Spacecraft Center 
NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 - CFSTI price $3.00 



PRECEDING PAGE BLANK NOT FILMED. 


FOREWORD 


The present report is one of a series of six reports, published simul- 
taneously, which describe analyses and computational procedures for: 1) pre- 

diction of the in-depth response of charring ablation materials, based on one- 
dimensional thermal streamtubes of arbitrary cross-section and considering 
general surface chemical and energy balances, and 2) nonsimilar solution of 
chemically reacting laminar boundary layers, with an approximate formulation 
for unequal diffusion and thermal diffusion coefficients for all species and 
with a general approach to the thermochemical solution of mixed equilibrium- 
nonequilibrium homogeneous or heterogeneous systems. Part I serves as a 
summary report and describes a procedure for coupling the charring ablator 
and boundary layer routines. The charring ablator procedure is described in 
Part II, whereas the fluid-mechanical aspects of the boundary layer and the 
boundary-layer solution procedure are treated in Part III. The approximation 
for multicomponent transport properties and the thermochemistry model are 
described in Parts IV and V, respectively. Finally, in Part VI an analysis 
is presented for the in-depth response of charring materials taking into ac- 
count char-density buildup near the surface due to coking reactions in depth. 


The titles in the series are: 


Part I 

Part II 

Part III 
Part IV 

Part V 

Part VI 


Summary Report: An Analysis of the Coupled Chemically Reacting 

Boundary Layer and Charring Ablator, by R. M. Kendall, E. P. 
Bartlett, R. A. Rindal, and C. B. Moyer. 

Finite Difference Solution for the In-depth Response of Charring 
Materials Considering Surface Chemical and Energy Balances, by 
C. B. Moyer and R. A. Rindal. 

Nonsimilar Solution of the Multicomponent Laminar Boundary Layer 
by an Integral Matrix Method, by E . P. Bartlett and R. M. Kendall. 

A Unified Approximation for Mixture Transport properties for Multi- 
component Boundary-Layer Applications, by E . P. Bartlett, R. M. 
Kendall, and R. A. Rindal. 

A General Approach to the Thermochemical Solution of Mixed Equilib- 
rium-Nonequilibrium, Homogeneous or Heterogeneous Systems, by 
R. M. Kendall. 

An Approach for Characterizing Charring Ablator Response with In- 
depth Coking Reactions, by R. A. Rindal. 


This effort was conducted for the Structures and Mechanics Division of 
the Manned Spacecraft Center, National Aeronautics and Space Administration 
under Contract No. NAS9-4599 to Vidya Division of Itek Corporation with Mr. 
Donald M. Curry and Mr. George Strouhal as the NASA Technical Monitors. The 
work was initiated by the present authors while at Vidya and was completed 
by Aerotherm Corporation under subcontract to Vidya (P.0. 8471 V9002) after 
Aerotherm purchased the physical assets of the Vidya Thermodynamics Depart- 
ment. Dr. Robert M. Kendall of Aerotherm was the Program Manager and Prin- 
cipal Investigator. 
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ABSTRACT 

A general equilibrium and nonequilibrium chemical state procedure is 
developed and applied mathematically to a number of open and closed thermo- 
dynamic systems. The conventional equilibrium relations are developed in 
terms of a set of base species. The base species concept is then extended 
in order to treat mixed equilibrium and nonequilibrium systems in a general 
fashion. The specification of controlling reactions is used to create non- 
redundant equation sets as equilibrium is approached. 

The treatment of open system mass balances within the basic framework 
of the state solution permits direct surface state calculations considering 
boundary- layer transfer relations and surface-condensed phase removal rela- 
tions. The relations defining the state downstream of an oblique shock are 
also included in the basic solution procedure permitting direct evaluation 
of these relations without recourse to subordinate iterative schemes. 

The factors affecting convergence within the framework of the Newton- 
Raphson iterative procedure are discussed and the techniques employed in the 
study are described. The means of evaluating state derivatives is described 
and relations presented for specific examples. 

The computer program which performs the equilibrium state solutions 
according to the methods presented, the Aerotherm Chemical Equilibrium (ACE) 
program is described briefly. The program which contains some of the nonequi- 
librium features of the analysis, the ACE/KINET program is also introduced. 

The current operational status of both programs is presented. 


v 
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A GENERAL APPROACH TO THE THERMOCHEMICAL SOLUTION 
OF MIXED EQUILIBRIUM-NONEQUILIBRIUM, HOMOGENEOUS 
OR HETEROGENEOUS SYSTEMS 

SECTION 1 
INTRODUCTION 

In the study of high-energy boundary-layer phenomena, thermochemical 
processes can be of dominant importance . This is particularly true when 
these boundary layers interact with chemically active surfaces. In the pres- 
ent study, interest is directed toward the superorbital reentry of the Apollo 
vehicle and the thermochemical response of its heat shield. The requirements 
for evaluating the chemical state of homogeneous and heterogeneous systems in 
this study are extensive. These requirements include the determination of 
the chemical state after normal or oblique shock wave compression, during the 
isentropic expansion of the inviscid shock-layer gases, within the boundary 
layer and at the chemically active surface. In the last two instances these 
state calculations are coupled with complex mass-balance relations. Many 
chemical-state solution procedures have been documented to treat reasonably 
standard closed systems, such as those associated with expansion processes. 

(See, for example. Refs. 1 through 4.) For open systems, only a few direct 
solution procedures have been documented. The surface state calculation for a 
general chemical system described in Ref. 5 falls into a special class of open 
system problems. Because of the number of requirements imposed upon the chemi- 
cal-state routines in the present study, the general treatment of a variety of 
chemical systems became a major effort. The inclusion of a general kinetic 
model, ionization, and the extensive bookkeeping associated with the downstream 
introduction of new species is of major importance in the formulation of the 
general problem necessary for thoroughly treating the coupled boundary-layer 
problem. 

Starting from reasonably fundamental relations, the procedures adopted as 
a part of this research effort are described in this report. The basic rela- 
tions are presented in the next section followed by sections on the solution 
procedure and the evaluation of thermodynamic properties. 

These techniques have been built in greater or lesser extent into the 
Aerotherm Chemical Equilibrium (ACE) program and certain special modifications 
of it. The final section of this report specifies the exact status of these 
routines and the extent to which the general formulation presented herein has 
been implemented. The procedures are presently limited to equilibrium, except 
that selected homogeneous reactions can be considered frozen within the bound- 
ary layer and selected heterogeneous and surface catalyzed kinetically con- 
trolled reactions can proceed at the material surface. 



SECTION 2 
BASIC RELATIONS 

In this section the relations which are required to specify the chemical 
state of a system are presented. Basically four types of relations can be 
considered in a general open system. These are the equilibrium relations 
applying to those reactions which can be considered as generally equilibrated, 
the nonequilibrium relations for those reactions which can be (but are not 
necessarily) out of equilibrium, the mass-balance relations, and those addi- 
tional state constraints imposed on the system. 

2.1 EQUILIBRIUM RELATIONS - TOTALLY EQUILIBRATED SYSTEMS 

In a chemical system there will exist, in the general case, a set of in- 
dependent equilibrium reactions. All other equilibrium reactions will be equiv- 
alent, both physically and mathematically, to this independent set. Consider 
for example the simple H, 0, HO, H^O, H 2# 0 2 , H 2 0 2 , O^ system. Six indepen- 
dent equilibrium reactions can be written, for example 


H + 0 


OH 

(i) 

2H + 0 

* 

H 2° 

(2) 

2H 


h 2 

(3) 

20 


°2 

(4) 

2H + 20 


H 2°2 

(5) 

30 

* 

°3 

(6) 

Other reactions such as 




H 2 + ho 2 


H 2° 

(7) 

are merely linear combinations of the six 
equations (i.e., Eq. (7) = Eq. (2) - Eq. 

(arbitrarily selected) 
(3) - ^Eq. (4)) . It can 

independent 
be shown 


that in a completely equilibrated system the number of independent equations 
is usually equal to the number of molecules less the number of elements.* The 


An exception occurs when two or more elements are in the same ratio in all 
molecules of a system, e.g., the system N0„, N^O. ^ as one ( n °t zero) indepen- 
dent equations . 
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modification of this relation for systems that are not completely equilibrated 
will be considered in Section 2.2. The selection of this set of indepen- 
dent reactions can be done arbitrarily, but it is convenient to establish 
some consistent technique. Most such techniques are based on the pre-selection 
of a set of species usually equal in number to the number of elements. The 
formation reactions of all other species from this base set represent the 
independent set of equilibrium reactions. The base species must be selected 
in such a fashion that no reaction can be written wherein reactants and pro- 
ducts are all base species. Thus in the above set HO and H 2 0 2 re P resent an 
invalid base set whereas HO and 0, HO and H, etc., represent valid sets. It 
has been reasonably common practice to select the monatomic gases as base spe- 
cies (Ref. 3, for example) since the formulation of the formation reactions 
is particularly convenient. There are advantages, however, in selecting a 
more general set, particularly when chemical kinetics are important. Consider- 
ing a set of base species N^, formation reactions for the remaining Nj species 
are of the form, indicated for example in Ref. 4, 



i 


where the are the stoichiometric coefficients of the formation reactions. 

In the preceding example with H and 0 as N-^ and N 2# respectively, the 
would be 


j 

1 

2 

3 

4 

5 

6 

7 

8 

i n. 

H 

0 

HO 

h 2 0 

H 2 

°2 

H 2°2 

°3 

1 H 

1 

0 

1 

2 

2 

0 

2 

0 

2 0 

0 

1 

1 

1 

0 

i 

2 

2 

3 


while with H 2 0 and 0 2 as and N 2 , the 


would be 


j 

1 

2 

3 

4 

5 

6 

7 

8 

i 

H 2° 

CM 

o 

H 

0 

HO 

H 2 

H 2°2 

°3 

1 h 2 0 

1 

0 

.5 

0 

.5 

1 

1 

0 

(N 

o 

CM 

0 

1 

-.25 

.5 

.25 

-.5 

.5 

1.5 
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(the atoms 


Mathematically, the v ^ ^ are obtained implicitly from the c ^ 
of element k in molecule j) by 


y c, . v . . 

L-> ki ji 


= c 


'kj 


(9) 


which for the latter table with j = 7 (i.e., H 2 0 2 ) 

for k representing hydrogen 2 x 1 + 0 x .5 

for k representing oxygen 1 x 1 + 2 x .5 

In matrix form 


indicates 
= 2 
= 2 



T 



C kj 


( 10 ) 


where the superscript T denotes transpose. Note that the square array 
c ki can be cons i dered as a subset of the larger rectangular array c^j . 

The set of independent formation reactions (Eq. (8)) for j ranging from 
+ 1 to N g , where N^ is the number of base species and N g is the total 
number of species, can be used to formulate a set of equilibrium constraints. 

At equilibrium, the second law requires that these independent reactions 
occur without change in free energy. Therefore 

a j - Iv G i (11) 

i 

where the Gj are the partial molar free energies of the species (also referred 
to as the chemical potentials) . The change in free energy is, by definition, 
equal to the isothermal reversible work performed by a steady flowing system 
in passing from one state to another. On this basis it is possible to relate 
the G j to the standard-state free energies, G? , that is, the free energy 
of the species at the same temperature but undiluted and at one atmosphere 
pressure. Thus 


G : - = 


J 


vdP 


( 12 ) 


-I isothermal, reversible 


where P° is one atmosphere. For a gas obeying the perfect gas law 


G? = - RT in p. 


(13) 
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if Pj is in units of atmospheres. Likewise, for an incompressible condensed 
phase containing only one species 




(14) 


When more than one condensed species coexists in a phase, the work of mixing 
must be included. For an ideal mixture* 



1 - p 

- RT in x 


(15) 


where 


is the mole fraction of species X in the mixture. 


In environments of general concern in high-temperature thermodynamics, 
Eqs. (13) and (14) are generally employed and, in addition, (1 - P)/p^ is 
assumed negligible in relation to the gas-phase work terms. On this basis, 
simplified equilibrium constant relations are obtained from Eqs. (11) and 
(1'3) for gas and condensed species (sub j or X) , 


Xn K 


AG? 


RT 


= x 


n Pj - E v ji * n p i 


(16) 


where the standard-state free energy change of the formation reaction for 
species j (or X) is defined by 



G° - V v . . G° 
3 La 3 ii 


(17) 


and the partial pressure of condensed species will be taken as one atmosphere 
in order that Eq. (13) will indicate no work of compression. 

The standard state free energy is a function of temperature only and is 
obtained for each molecular species from 


G° = H° - TS° 


(18) 


where enthalpies are obtained relative to some chemical base state, often the 
elements in their most natural form at 298°K and one atmosphere (JANAF base 


The ideal mixture assumption is satisfied by a mixture above which vapor pres 
sures are proportional to mixture mole fractions and whose vapors obey the 
perfect gas law. 
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state) , If any other base state is consistently adopted, the AG? will be 
unaffected. 

The stationary condition of the free energy at equilibrium expressed in 
Eq. (11) is consistent with the minimum free energy statement often utilized 
in seeking the equilibrium state. Although the formulation followed here 
differs from those followed in free-energy minimization approaches, the ulti- 
mate numerics can reduce to an identical iterative solution procedure. 

The solution of the set of algebraic equations (Eq. (16)) must be con- 
sidered in conjunction with other constraints including the pressure balance 

- * < i9 > 

j 

where the summation is over all gas-phase species. The detailed solution pro- 
cedure will be considered only after all relations have been presented. 

2.2 MIXED EQUILIBRIUM-NONEQUILIBRIUM RELATIONS 

When some reactions fail to equilibrate it is necessary to approach the 
selection of the independent set of equilibrium reactions with greater caution. 
In a general chemical system, certain sets of molecules can be treated as al- 
ways equilibrated. Between these sets certain independent equilibrium and 
kinetically-controlled interchange reactions may exist. A procedure for treat- 
ing mixed equilibrium-nonequilibrium systems is presented in this section. 

The following rules are established in order to organize the logic: 

1. Every species is assigned to one and only one set of molecules. 

2. A set may contain as few as one species. 

3. Each set has its own base species, i.e., that minimum number of 
species from which all other members of the set may be formed. 

4. Within each set all possible reactions between member species 
are equilibrated. 

5. Equilibrium interchange reactions involve species from more 
than one set. 

Consider the eight species of the 0-H system considered in Section 2.1? 

0, H 2 0; H; H 2 ? 0 2 , 0 3 ? HO, H 2 0 2 where five sets are divided by semicolons. 

For these sets, the following base species are appropriate: 0, H 2 0? H; H 2 ? 

0 2 ? HO where only the first set requires more than one base species. At 
this juncture only two independent equilibrium reactions have been formulated, 
namely 
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( 20 ) 


I - 5 ° 2 * ° 3 

2 HO ^ H 2 0 2 (21) 

Two independent equilibrium interchange reactions might be included in this 
system, for example 


H + OH H 2 + 0 

(22) 

0 + OH Oj + H 

(23) 


The effect of these reactions is to reduce the number of base species by two. 
For example H 2 and 0 2 can be deleted. The remaining base species and the 
array of formation reactions coefficients, are therefore 


\ 

3 

1 

2 

3 

4 

— 

5 

6 

7 

8 

i 

s \ 

0 

h 2 0 

H 

OH 

H 2 

°2 

H 2°2 

°3 

1 

0 

1 

0 

0 

0 

-1 

1 

0 

1.5 

2 

H 2° 

0 

1 

0 

0 

0 

0 

0 

0 

3 

H 

0 

0 

1 

0 

1 

-1 

0 

-1.5 

4 

OH 

0 

0 

0 

1 

1 

1 

2 

1.5 


in the general relation 


Z'-jA = N j < 8 > 

i 

To illustrate, the equilibrium formation reaction for 0 3 (i.e., j = 8 in 

the above table) from the base species is given by 

3/2 0 + 3/2 OH - 3/2 H & 0 3 (24) 

If it is assumed that all other interchange reactions are frozen, the 
formulation of the equilibrium-nonequilibrium aspects of the program are com- 
plete. In the totally equilibrated chemical system, conservational constraints 
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are often applied to the elements. in the system just presented, however, 
additional conservational constraints are required. In general these con- 
straints take the form 


,25 > 

j 

where is the number of moles of species j in a unit mass of system and 

cu is a conserved variable relating to the "elemental" composition of a unit 
mass of the system. One such constraint is applied for each base species. To 
check the validity of this relation, consider the effect of a unit reduction 
in n. , . According to Eq. (8) the n^ will increase by v j • ^ • Thus in 

Eq. (25) two terms are modified, namely, v^n^ and v j'i n j l * s ^ nce = 

1.0, it is apparent that the increase in the latter term exactly balances the 
decrease in the former. In effect, the base species become the elements of 
the system, and their total masses can be treated as the conserved variable of 
the system. This generalized concept of the conserved "elements" of the sys- 
tem is extremely important to the present development and thus further examples 
are appropriate. In the completely equilibrated system previously presented, 
it is easy to accept that 


j n j =a k (26) 

j 

where , the total number of atoms of element k in a unit mass, is a con- 

served variable. If we premultiply both sides of this expression by the in- 
verse | . J ~ ^ ' there results from Eq. (10) 

Z v ji n i = hkil" 1 a k = a i (27) 

i 

Therefore if the a^. are conserved then the will be conserved. In the 

system of Eqs. (1) through (6), this is tantamount to naming a new set of 
elements which are, in effect, the base species. The term "element" (in quotes) 
will be used henceforth in this report to refer to those atoms or groupings of 
atoms (i.e. grouped according to the base species formulae) which according to 
the equilibrium relations are conserved. 

In the general nonequilibrium system, certain kinetically-controlled 
reactions will be important. For example, in the H-0 system 


H + H + M-+H 2 + M (28) 

0 + 0 + M— 0 2 +M (29) 

H + H + 0 — ► H 2 0 (30) 
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are three reactions of possible interest, M being any third body. The rates 
of these reactions can be related to the partial pressures of the reactants and 
products, the ultimate equilibrium constraint appropriate to the reaction, and 
the kinetic coefficient. With the general kinetic reaction in the form 


z ^ - z ^ 


(31) 


its rate can be expressed generally by 


R = k_ 
m F 

m 


TT - - 


j K 


TT ph 3 


p 

jm 


L 3 


■m 3 


(32) 


or by 


R = K 


m 


£ 


exp ) l n p - exp 


£* - 


p . - in K 

3 P, 


(33) 


m 


where 


in K 


m 



RT RT 


z 



R 

jm 


G° 


(34) 


The net effect of these reactions is the modification of the "elemental" makeup 
of the system. The kinetic reactions will cause a net increase rate (moles per 
unit volume) 


r iT ^ ( ^jm ' ^jm )v ji R m (35) 

of "element" i . It is this relation which will be introduced into the 
conservational equations in order to establish the local state of the reacting 
chemical system. Certain additional data must be provided in order to perform 
the kinetic calculations. First, the specification of the stoichiometry serves 
to establish not only the effect of the reaction but also the reaction order. 
For example 


c* + 0 2 — CO 


(36) 
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and 


2C* + 0 2 — 2C0 


(37) 


are equivalent stoichiometric ally, the asterisk designating condensed-phase 

material. However, in Eq. (32) , the former relation results in a half order 

reaction (p + = 1.0)*, vdiereas the latter results in a first-order reaction. 
C 

The forward rate constant k F , hopefully based on experimental data, is 
represented with an Arrhenius m type function 



B exp 
m * 


E /RT 
a 

m 


(38) 


where the exponential factor establishes the probability of a collision having 
energy in excess of the activation energy and the factor B m represents a 

multitude of phenomena associated with the probability of success of a single 
collision (e.g. collision orientation). 

When kinetically-controlled reactions approach equilibrium, difficulty is 
often encountered in the treatment of the relevant conservational equations. 

To understand the nature of this difficulty, and thus the means of avoiding 
it, it is instructive to consider the simple time dependent character of the 
H-0 system previously described. Recalling that cu represents the moles of 
"element" i in a unit mass, and that r^ represents the rate of production 
of moles of "element" i per unit volume, it follows that 


dcu r^RT 

dtT = "W 


(39) 


At this point it is necessary to introduce a new concept. From the base 
species, a subset of base-base species can be obtained much as if all speci- 

fied kinetic reactions were permitted to equilibrate. For this example, 0 and H 
will be selected for this honor and the "formation reaction" for the remaining 
base species written as 


2H + 0 — ► H 2 0 

(40) 

H + 0 — ► OH 

(41) 


* ... 

As in the equilibrium relation, it is convenient to consider the partial pres- 
sures of condensed species as one atmosphere as a device to include hetero- 
geneous reactions in the same framework as homogeneous reactions. 
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These reactions will equilibrate if the kinetic reaction (Eq. (30)) and the 
reaction of either Eq. (28) or Eq. (29) have infinite rates.* More generally, 
Eqs. (40) and (41) can be written as 


I 




— N. 


l 


(42) 


where the index k represents the base-base subset of base species. it can 
be shown that 


L 


ik 


r . 

i 


= 0 


(43) 


which implies that atoms and/or certain combinations of atoms cannot be created 
or destroyed by chemical reactions in this system. Since there are two base- 
base species, two of the Eqs. (39) can be replaced by 


Defining 



1 


it follows that 


(44) 


(45) 



(46) 


For each base species i which is also a base-base species k , an 
equation of the form of Eq. (46) will replace the corresponding one of the 
form of Eq. (39) . The other Eqs. (39) are maintained unaltered in the system 
of equations and still contain the kinetic expressions. These will be referred 
to as the reactive mass balance equations. In the general case, a given 
reaction, m , will affect more than one of these equations introducing a term 
of the form 


The number of base-base species can still exceed the number of atomic elements 
in the system if insufficient kinetic paths have been specified to permit full 
system equilibrium. 
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( 48 ) 


If the first term in braces approaches infinity, the balancing of the overall 
reactive mass balance equation will reduce to setting the second term in 
braces equal to unity, or equivalently, 

" £nK P m = 0 (49) 


which is simply the appropriate equilibrium relation. If this reduction for 
reaction m occurs in more than one mass balance equation, a net loss of non- 
redundant relations can occur.* To avoid this, it is important that the reac- 
tive mass balance equations be combined in such a manner that the production 
terms from each near-equilibrated reaction is entered into only one equation. 

The means of establishing this combination of the reactive mass balance equa- 
tions is based on the selection of controlling reactions equal in number to 
the number of reactive mass balances. (This number equals the difference between 
the number of base species and the number of base-base species.) in turn the 
selection of the controlling reactions is based on the array Q^ m where 


^im 


I 


(p? - p R ) v . ■ 


exp 


V F 


R * 

jm* n Pj 


(50) 


The performance of a conventional pivotal Gaussian reduction on this rectangular 
array results in the selection of certain of the Q ^ as pivotal terms. it is 
the column numbers, m , of these terms which are considered as the indices of 
the controlling reactions, one having been selected for each reactive mass 
balance, i . The combination of the reactive mass balance equations such as 
to eliminate terms consequent to the controlling reactions from all but one 
equation thus becomes a relatively straight-forward manipulation following the 
same steps as the original Gaussian reduction. This manipulation is carried 
forward independently for both the mass balances and the kinetic relations, 
before the two sets are combined, in order to avoid the loss of important 
significant figures. In actual practice, this manipulation is performed in 
conjunction with the transformation indicated by Eq. (44) and is effected by 
an augmented transformation matrix. 

*Such a loss is only a result of the limited significant figures retained in the 
calculational process since differences in the other terms in the equation 
could still be used if unlimited significant figures existed. 
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2.3 MASS BALANCE RELATIONS 

In the preceding sections the equilibrium and nonequilibrium relations 
have been developed for a fairly general chemical state. These relations are 
in themselves insufficient until other relations, in particular the mass 
balance relations, are imposed. In the case of kinetic control, the time 
dependence of the system must be equated to flow rates and other rate dependent 
parameters entering the mass balances. Likewise, in diffusional systems the 
local state is determined by mass-balance relations associated with mass- 
transfer processes. in the following subsections the mass balance relations 
appropriate to various systems will be presented. 

2.3.1 Expansion of Isolated Systems 

In the expansion of a fixed mass, closed, adiabatic system it is usually 
appropriate to trace its state history as a function of static pressure. If 
the process is reversible, the entropy is constant and the local state is not 
a function of the time history (process path) of the expansion. Such systems 
satisfy the simple mass balance constraint 

ou = constant (51) 

This equation implies either total equilibrium or a mixed equilibrium- frozen 
chemical process. If, however, finite reaction rates are important, the path 
ceases to be reversible, entropy rises, and the time history of the expansion 
must be considered. If the pressure is a known function of time, the expansion 
can be treated as 


state = f(au, s * p ) 


(52) 


where the state includes such terms as dcu/dt and ds/dt . The rate of change 
of the "elemental" composition is obtained from Eqs . (33) and (35), whereas the 

rate of increase in entropy (see, e.g.. Ref. 6, Eq. (3.32)) is obtained from 


ds V V- RTx 

dt T l m P ?n 

m 


(53) 


Where A^ (the "affinity") is defined as 


A m = RT 



I 



djm )£n Pj_ 


(54) 
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Thus 


ds 

dt 


- 


L m 


(}i P - p R ) £n p , 
J m ]m 


R El 

m 


(55) 


This derivative is well behaved, even as equilibrium is approached, and may be 
evaluated explicitly if desired. The derivatives, however, must be 

treated implicitly if any hope for near equilibrium solutions is to be main- 
tained. Once a particular formulation is adopted, the techniques suggested in 
the nonequilibrium presentation (Section 2.2) can be introduced in order to 
assure consistent solution validity. Because of the simplicity of the mass- 
balance relation for this process, it is practical to include the kinetic mass 
balances directly with the iterative solution of the chemical state. This 
state calculation includes the relations previously presented together with the 
entropy constraint, namely 

£ Pj S. = P^s (56) 


where the entropy of a perfect gas species j is related to the standard 
state entropy (at one atmosphere pressure) by 


Sj = S° - R in Pj (57) 

and S° is a function of temperature only. For condensed species the simpli- 

3 Q 

fying assumption presented previously leads to . In all standard 

conservation equations, the partial pressure, p^ , assumes the general defini- 
tion of the number of moles of species j for P gas -phase moles. 

2.3.2 State Calculations for Open Systems 

The evaluation of the state in a general open system involving diffusive 
and convective mass and energy fluxes is most generally performed as a sub- 
ordinate solution. For example, in the boundary-layer solutions of current 
interest, state solutions are required at several interacting locations. Thus, 
in Ref. 7 (Part III of the present series of reports) state solutions are 
required based on assigned "elemental" mass fractions (or cu ) , enthalpy and 
pressure, i.e. 


state = f(a i , h, P) 


(58) 
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This solution provides to the boundary-layer solution the detailed state 
including thermodynamic, transport and radiative properties as well as the 
production rates dcu/dt . The last term must be included with the general 
mass-balance relations of the boundary-layer program. 

The specific relations used to achieve the state solution are the equi- 
librium equations (Eq. (16)), the mass balance relations (Eq. (25)), the pres- 
sure constraint (Eq. (19)) and an enthalpy constraint 

^ Pj H. = P% (59) 

j 

which involve no greater complexity than the conventional isolated system 
equilibrium solution. The coupling between this solution and the boundary- 
layer solution requires not only the evaluation of the production rates, 
da^/dt but also the rate of change of these rates with respect to the 
independent parameters on the right- hand- side of Eq. (58) . The rates are 
determined with Eqs. (35) and the derivatives by use of relations to be 
developed in section 3 of this report. 

Again, the problem of stability threatens when equilibrium approaches. 

However, by following the approach previously presented, the kinetic terms can 
be treated separately while all other terms of the mass-balance equations are 
being collected. These equations can then be rearranged and combined with the 
kinetic relations in such a way that controlling reactions again affect only 
one equation at each location. Because of the overall implicit character of 
the boundary- layer solution, this procedure will, on convergence, yield valid 
stable solutions, it has been found, however, that the introduction of equi- 
librium type relations into the set of boundary- layer equations can destroy 
the linearity of the system. Therefore the approach of equilibrium by the 
kinetic equations included in the boundary- layer mass balances will probably 
delay convergence and necessitate the inclusion of certain iterations con- 
straints . 

2.3.3 Surface State Solutions 

A more complex set of mass-balance relations are introduced when surface 
state solutions are sought. Coupling between boundary layer, internal conduction, 
and surface mechanical removal solutions may be involved in these relations. in 
effect, all the other mass-balance solutions become subordinate to this solution. 

Two types of boundary-layer representations are considered in this section. 

One is the transfer-coefficient correlation of mass transfer using the Z*-potential 
developed in Ref. 5 and also in Appendix B of Ref. 8. In Ref. 7 the fluxes at 
the wall are expressed using an influence coefficient approach. Mathematically 
these approaches reduce to 


h = p e u e C M^ ‘ 


(60) 
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and 


= 




(h 


h°) 


9j. 


3 (pv), 


(p v) 
r w 


(pv) 


w 


kdK k 




( 61 ) 


respectively. The nomenclature is as given in those references and the list 
of symbols of this document. The diffusional mass flux of "element" i 
is coupled with the surface mass-balance equation 


m K 
c c . 


iti K 

g g* 




m K, . - (pv) K. 
r^ U v K w i 


= 0 


(62) 


This equation includes fluxes due to the pyrolysis gas generation rate, m , 

y 

the char consumption rate, m c , the surface mechanical removal rates of spe- 
cies i due to various mechanisms, m r ^ , the bulk gas-phase convective flux, 
(pv) , each times their respective fraction of "element" i , and the diffu- 

sive flux, j., of "element" i . An overall mass balance can also be written 

i _ 

by summing Eq. (62) over all i and noting the 2_,j^ = 0 , namely, 

® c + ® a - Z"* r “ (P v ) w = 0 ( 63 > 

C y l l 


Because of the different treatment used with Eqs. (60) and (61) # they will 
be developed individually, but generalized in a consistent format. 


2. 3. 3.1 Surface Transfer Coefficients 

In many boundary- layer analyses, the net diffusional flux of mass, heat 
and/or momentum to or from the wall are the only results of practical concern. 
Because of this, transfer coefficients have become extremely effective tools of 
analysis when it is possible to develop good correlations by which their values 
can be adequately estimated. Using the transfer potentials developed in Ref. 5 
and Appendix B of Ref. 8, a directly coupled equilibrium surface thermochemistry 
routine has been developed. Combining Eqs. (60) and (62) and rearranging yields 


1 
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p e u e C M 


K 


HI 


p e U e C M 
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= a 


z. 

1 




<P v >w ~ V % 

„ „ - V- K, +L— 7-S 
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H e e M 
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M e e M 
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■ii 


(64) 
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where the JL subscript refers to the possible condensed species in the system. 
The mass fraction of "element: i in the species i, K^, can be expressed as 


7I\. 

*Xi = V 4i ^ 


(65) 


Similarly, the gas phase mass fraction of "element" i can be expressed in 
terms of the partial pressures of the gaseous species, j, namely, 


*i P»L ZZ V ji P j\ 


( 66 ) 


where W{ is the mean molecular weight of the gases adjacent to the wall, 
that is. 


PAL 


■ Z 


(67) 


where the summation over j is for gaseous species only. Recall that the 
are the stoiciometric coefficients in the formation of species j from 
"element" i or the atoms of "element" i in a molecule of species j. 

It is convenient to define the term p^ for condensed species in a spe- 
cial fashion for this open system mass balance, namely 


P ‘ »eWj! 


(68) 


Combining Eqs. (63) , (65) , (66) and (68) with Eq. (64) yields 


Z f + lit 

a = _i + _L_ -2 £ 

Z. 71 \ . P57? p u C M 

i i g \ r e e M 


E p i 7n A 


VIA 


g I 


z 




V ji P j + 


(69) 


From Ref. 5, the definition of z^ is found to be 


s . . ^y P .m 

i pm 2 L-j Fj 


(70) 
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where F* is equal to F ^ raised to some power, usually 2/3 for boundary 
layer applications. The diffusion factor F ^ is described in detail in 
Ref. 5 and in Ref. 9. Also 




Y 

p j PS 


(71) 


Combining this result with Eq. (69) yields, after introducing a normalizing 
parameter F (= 




\ + “ c Ev/V V V 

fc: p kt\L Pj^ji + p i v ^i + L 


v ii 

p. -Ui- (72) 
. ^ F*/F 


From Eqs. (67) and (71) , F is defined by 


j , i , ^ 

EPj^jA^ 


(73) 


One of the more elusive aspects of surface-state solutions is the ade- 
quate specification of the mechanical-chemical surface constraint. The pres- 
ent formulation is based on the following set of contraints for condensed- 
phase species. 


p* 


0 if T < T 

F l 


(74) 


and 


E v ii * n p i * in K Pje (?5 > 

i 


with the equality applying to one species with T T . The first equation 

*1 

implies that a particular condensed species cannot leave the surface until the 
surface temperature is at or above that species fail or flow temperature, T_ 

The second relation states that all present condensed species are in equilibrium 
with the base species. The inequality applies to non-present condensed species 
and prohibits a super-saturated vapor state. This is equivalent to saying 
that at 100°C 
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(76) 


p H 2 0 * 1 atm. 


The requirement that one species be at or below its fail temperature establishes 
the structural limitation of the surface. A typical result might show a surface 
at 2500°K with the equilibrium holding for SiC*and SiO*, but if SiO* has been 
assigned a fail temperature of, say, 2300°K, P siQ and thus m siQ will be 


positive indicating liquid removal of SiO^ . Ihe SiC*, with a fail temperature 
greater than 2500°K, represents the surface constraining species. 


2.3.3. 2 Boundary Layer Coupled Surface State 

In order to couple the boundary-layer solution to the surface state equa- 
tion, several approaches are feasible. Because of the nonlinearity of surface- 
state solutions, there are significant advantages to isolating this solution 
from the more linearly behaving transfer equations of the boundary-layer. The 
chemical -state routines expect and can accept quite nonlinear systems of equa- 
tions. The boundary-layer aspects of the surface- state calculations are thus 
reduced to influence coefficients, similar in form to the transfer coefficients 
discussed in the previous subsection, which are introduced into the surface 
state equations. This linear reduction of the boundary-layer equation is 
discussed in Ref. 7. With this approach in mind, the second surface diffu- 
sional flux equation (Eq. (61)) can be introduced into Eq. (62) yielding, 
after some rearrangement 


% 


rn.K_ + mJL + 4^t-(pv)° + h° K? - j° 

ah k 3K^ K 1 


c c. 


i g g ± a(pv) w ^ v 'w 


a i “ \ 


5jj 

3 (pv) w w 1 dh 


3j i V 3j i ~ , „ ~ V 

+ TvT* + K. + (Pv) K. + L 1 


K. + (Pv) K. + L™ r 
k 31^ k w 1 t *1 £1 


(77) 


where the partial derivatives and prior values (°) of the variables are 
developed by the boundary-layer solution procedure. 

Making use of Eqs . (63), (65) and (66) and introducing the following 

definitions 
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(78) 
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+ y h 1 I P 3 H i + ^ u ii p '« 


(82) 


Introducing a new coefficient 


... = y 

J 1 L—> 




(83) 


reduces Eq. (82) to 


p Vi = 


m + m 




p Vf. + 

^ l 



* Z p i x ji * y h Z p j“j 

j 


(84) 


This equation is similar to Eq. (72) and can be reduced to that equation if 


m + m 

m + m is replaced by —2 — — - 

9 c p e u e C M 


(85) 


Y-l and are replaced by 0 


(86) 


v . . F 

is replaced by — — 


(87) 
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The inclusion of heterogeneous kinetics in these equations is accomplished 
by adding the production rate of "element" i per unit surface area to the 
left side of Eq. (72) or (84) . The subsequent treatment of this tern is quite 
similar to that previously presented for homogeneous kinetics. Thus with the 
k of Eq. (38) defined in terms of stoichiometric moles per unit of time and 
surface area, the production term to be appended to the left of the final 
Eq. (84) is with r^ defined according to Eq. (35) . The same term 

should be added to Eq. (7 2) but the k p are then normalized by p u C . 

As in the preceding kinetic problems, a set of controlling reactions is de- 
termined and the mass balance equations are rearranged so that each of these 
reactions appears in only one of these equations, at which juncture the kinetic 
terms are introduced into the mass-balance equations. 

2.4 OBLIQUE SHOCK RELATIONS AND SUMMARY OF STATE RELATIONS 

Certain constraints in addition to mass and equilibrium balances have 
been mentioned in the preceding paragraphs, it is well to summarize these 
briefly and to introduce one additional set of constraints, namely those 
associated with flow across a shock wave. Up to this point, total pressure 
has always been an independently assigned variable, thus 



j 


is an obvious constraint on the partial pressures. Other state constraints 
have also been mentioned. In addition to temperature specification they are 


and 


£ PjHj = hP^ (89) 
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(90) 


for specified enthalpy, h, and entropy, s, both defined per unit mass. 

In the case of an oblique shock wave, none of the above constraints apply. 
Instead, conservation of energy, mass and momentum are required. Using the 
subscripts, 1 and 2, to denote conditions upstream and downstream of the shock, 
respectively, with 0 as the angle between the flow vector and a normal to 
the shock wave, these equations are 


energy 




( 91 ) 
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mass 


P 1 u 1 cos0 1 = P 2 U 2 coS ^2 


normal momentum 


Pl + p 1 u 1 cos0 1 u 1 cos0 1 = P 2 + p 2 U 2 COS0 2 U 2 COS0 2 

tangential momentum 

p ^u^cos0^u^sin0^ = p 2 u 2 cos6? 2 u 2 sin0 2 

Combining Eqs . (92) and (94) yields 


u^sinG^ = u 2 sin0 2 


that is, the tangential component of velocity is preserved across the 
Eqs. (92) and (93) combine to yield 


Pi + 


(p 1 u 1 cos0 1 ) 


(p 1 u 1 cos0 1 )' 


i ‘ ^ + e 2 


Now 


u 2 “ ( u 2 sin0 2^ 5 + ( U 2 cos ^ 2 ^ * 


which with Eqs. (92) and (95) becomes 


U 2 = ( u 1 s in0 1 ) 2 + 


p 2. u i cos ^i ' 


With this, relation (95) becomes 


h l + 2 


— = h 2 + j ± (p^cose.^ s 
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shock. 

(96) 

(97) 
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Presuming knowledge of the upstream conditions, Eqs . (96) and (99) contain 

as unknowns only and These two equations can be further re- 

duced to 


Z p j + (p i u i cos( i> 3 H 


p i + 


(P 1 u 1 cos0 1 ) i 

Pi 


Z p j H j + 1 (P 1 u 1 cos0 1 ) s 


(RT) 3 


P % 


h l + 


(UjCOsQp 3 
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( 100 ) 


( 101 ) 


where the non-subscripted variables are downstream of the shock. The first 
of this pair of equations replaces the more conventional pressure constraint 
and the latter the enthalpy constraint. 

2 . 5 SUMMARY 

In this section an attempt has been made to formalize the basic relations 
so as to simplify the generation of an orderly solution. Unfortunately dealing 
with nonlinear equations such as these is never straightforward and is subject 
to many pitfalls. In the next section, the procedures as adopted in the cur- 
rent solution technique are described. 


SECTION 3 

SOLUTION PROCEDURES 


3.1 INTRODUCTION 

The solution to a set of simultaneous nonlinear algebraic equations can 
be either trivially simple or agonizingly difficult, depending on the lin- 
earity of the system and the depth of coupling existing between the equations . 
None of the problems formulated in the last section fall into the first class 
and some fall into the latter. The basic formulation adopted is relatively 
conventional, and will be described first, followed by some discussion of the 
pitfalls that can be encountered and devices adopted to circumvent them. 

3.2 BASIC FORMULATION 

The most direct method of solving a set of nonlinear algebraic equations 
is the Newton-Raphson procedure. Its application is straightforward in con- 
cept but in reality many choices occur during the formulation of a specific 
problem, choices which can affect the success or failure of a specific solu- 
tion. The method itself is the extension of Newton's iterative method to 
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multi-dimensional problems. Errors are evaluated for each of the equations 
based on a set of trial values for the unknown independent variables. The 
rates of change of these errors with respect to these independent variables 
are analytically determined, also based on the trial values. The differential 
relation 


dE, dE, 

dE l - dV l + aVj ^ + ••• (102> 

will apply for each error and the set of independent variables V^, ••• • 

Considering this as a set of simultaneous equations in dV-^, dV^, etc,, solu- 
tion is readily obtained as 


dV, dV ? 

dv i = 917 “1 + sSJ dE ; 


(103) 


where the partial derivatives are the elements of the inverse of the matrix 
of the partial derivative coefficients of Eq. (102) . Presuming these deriva- 
tives to be constants, the necessary corrections to each independent variable 
could be obtained as, for example 


dV, dV ? 

AV i = AE i + AE : 


(104) 


where the AE^, AE 2 , ... are simply -E^, _E 2 ' ••• if the errors are to be 

driven to zero. 


It should be noted that if some function of V is substituted for V 
in the above equation 


Af (V 1 ) 


df (vp 


dV, 


3V. 3V. 

3E^ AE 1 + 3E^ AE 2 + 


(105) 


implying simply that in the formulation of the derivatives no commitment is 

made with regard to the optimum means of expressing the corrections of the 

independent variables, for example, in terms of in p^ or simply p ^ . In 

the formulation followed here, in p., p. , in T, and in(P^) are taken 

J *• 

as the set of independent variables, but corrections are often in terms of 
Pj, 1/T and P^. In some systems this choice yields linear mass balance 
equations which if once satisfied will never deviate. 
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In Table I the relations developed in the preceding section are summarized . 
In addition to the relations, given as error equations, the derivatives of the 
errors with respect to each of the set of independent variables are given. It 
should be apparent that many arbitrary choices go into the relations for these 
derivatives. As an example of this arbitrariness, consider the closed system 
mass balance and the derivative with respect to in(P^) . In the table 
is given but the actual program permits either that term or the summation of 
v jiPj over all j and i. At convergence these terms are identical, but 
during the convergence process different paths are followed by the two expres- 
sions. Cases have been encountered where convergence was unsuccessful with 
the first expression and successful with the latter. 

3 . 3 SOLUTION CONVERGENCE 

In general the convergence of the set of equations appropriate to a 
particular problem will depend on a number of factors in addition to the for- 
mulation of the derivatives. Correction coordinates, initial estimates, and 
correction restraints are all major factors. The following subsections will 
explore these factors in some detail. 

3.3.1 Correction Coordinates 

It has already been mentioned that the choice of coordinates can demon- 
strably effect the convergence of a system. Obviously a coordinate system 
that most nearly linearizes the relations between errors and independent 
variables is to be sought. Thus we note that with the closed system mass 
balance relation, the error expression is linearly dependent on p %, p ^ , and 
p^ 9 thus making them prime candidates as the independent correction variables. 
Likewise, the gas-phase equilibrium equation with error expressed as shown is 
linear in 1 /t and in p^. , creating a bit of conflict between the two expres- 
sions for p j . 

Another specific example involves sets of species which are significant 
in only one mass balance, in particular ions which may be consequential in 
only the charge neutrality balance. For singly charged or ionized species, 
it can be shown that both charge balances and equilibrium can be expressed 
linearly in terms of the product PjPj where I denotes the base species 
associated with charge neutrality (usually the electron) and j represents 
all charged species. The mass balance (i.e., charge neutrality balance) is 
simply multiplied by p . The equilibrium relations (Eq. (16)) are written 
as 
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where the product is for all the neutral base species. For singly charged 

species with the electron as I, Vj I is either +1 or -1 and the linearity 

of p p . is obtained. (Note that p?. is included in p p..) This repre- 
x j x I 3 

sents a considerable departure from the previous evaluation of equilibrium 
errors, this error equation being related to the previous one (Table 1) by 

Ej =-p I Pj £exp (Ej) - lj (107) 


The purpose of this discussion has been to establish the flexibility of the 
Newton-Raphson procedure with regard to dependent and independent variable 
specifications and the advantages that can be achieved by seeking a nearly 
linear dependence of the errors on the independent variables. 

3.3.2 Correction Restraints 

In this highly nonlinear application of the Newton-Raphson technique a 
variety of constraints with regard to independent variable corrections are 
necessary. These constraints all manifest themselves in a damping factor 
which limits the extent which the solution is advanced down the correction 
vector. 

This factor is determined by use of the following limiting relation, 

4 + 4 (£n Pj) 

5)>hi ‘ < 108 > 

3 - Un p.) 

J io 

that is, if a predicted increase in p. from a value p_: yields a value 

_ 3 J .1 o 

of p-; exceeding that indicated by the right hand side of Eq. (108) , a 
J hi 

damping factor will be applied such as to achieve the equality. For the 
case of decreasing p ^ this equation can be manipulated to 

-4 + 3 {in Pj) 

Un p.) ^ = 1 (109) 

J Jto 4 + (in p .) 

3 hi 

with the same implications. The following table demonstrates the nature of 
the equality in these relations. 


26 


in Pj 

JLo 

— CO 

-13 

-5 

-3 

-2.33 

-2 

-1 

0 

hi 

-4 

-3 

-2 

-1.33* 

-1 

-.6 

0 

1.33 


The term, p ^ , is the partial pressure of species i normalized by its 
relative contribution to the mass balance equation wherein it is most signi- 
ficant. Specifically, 



max over 
all j 


max over 
all k 


(HO) 


where v i s t ^ ie sum of a ll terms involving p. in mass balance k. . -yr 

J / -K 3 

gas phase species whose initial value of p^ is less than the appropriate 

starred entry in the above table a logarythmic partial pressure correction is 

applied. Otherwise linear corrections are employed. After evaluating all 

corrections the minimum damping factor is determined, and is subsequently 

applied uniformly to all corrections. 


3.3.3 Initial Guesses 

It is obvious that a good first guess can save time in any iterative solu- 
tion. In the present formulation these guesses are generally based on previous 
solutions and only the initial stagnation or shock solution does not have the 
benefit of prior solution. This solution is readily obtained from practically 
any first guess, since the stagnation state is usually at relatively elevated 
temperatures and has a fixed "elemental" composition. In the subroutine ver- 
sion of the chemical-state program used in conjunction with the boundary-layer 
procedure, first guesses are generally based on solutions at the same boundary- 
layer transverse location stored during prior iterations in the boundary-layer 
program or from solutions at the preceding axial station. 

Because of the introduction of new species by the wall material it is 
necessary to initialize their compositions when the corresponding elements 
appear in the state solutions. Likewise if a species disappears, e.g., as 
the edge of the boundary layer is approached or because the sequence of bound- 
ary layer iterations results in the termination of surface mass addition, it 
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is necessary to zero the species in a fashion that will not result in a singu- 
lar solution for the rest of the equations. 

It is apparent that bookkeeping becomes a major factor in the state 
programs if efficient and stable repetitive utilization is to be made of the 
routines. This bookkeeping establishes optimum first guesses, determines 
which atomic elements are present, and zeros or initializes the appropriate 
molecular species . 

3.4 SOLUTION MECHANICS 

Up to this point the more basic aspects of the solution procedure have 
been presented. The actual mechanics also need to be described. The computa- 
tional procedure is composed of eight major parts, the bookkeeping, the input 
and data organization, the evaluation of temperature-dependent thermodynamic 
properties, the evaluation of errors and the derivatives of the errors with 
respect to the independent variables, evaluation and integration of kinetic 
terms into the arrays of errors and derivatives, the inversion of the reduced 
forn of the resultant matrix, the evaluation of corrections as limited by the 
constraints, and the determination of properties and property derivatives 
(after final convergence) . 

These eight parts are represented by eight routines in the Aerotherm 
Chemical Equilibrium (ACE) program. Each of these will be briefly described 
in the following paragraphs. 

EQUIL 

This routine is the main routine of the ACE program (or subprogram in con- 
junction with the CABLE program) . It controls the majority of the bookkeeping, 
develops the constant terms of the mass balances, outputs the solution, con- 
trols the main iteration, and exercises certain limited solution constraints. 

INPUT 

This routine controls the reading of element and species data, the selec- 
tion and setting-up of base species, the evaluation of the array, flag- 

ging of condensed species, and the setting-up of the very first guesses. 

THERM 

The evaluation of such temperature-dependent thermodynamic properties 
as enthalpy, entropy, specific heat and free energy are determined for each 
species by this routine. 
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MATER 


For the specific option being considered, this routine evaluates all 
errors and error derivatives according to the relations given in Table I . 

It also prepares a matrix of coefficients based on a reduced set of indepen- 
dent variables. This is accomplished by substituting the expression for 
Xn Pj obtained from the equilibrium relations (Eq. (16)) into the mass 
balance relations. This simple algebraic substitution produces a very sig- 
nificant reduction in the order of the coefficient array. Thus in the mass 
balance equations 


dE 

a in Pj ' d £n Pj (HD 


is replaced by 


dE. 


d Xn p . 


t 

L i 


v . . d Xn p. 
Di 



Xn T 


- dEj 


( 112 ) 


The variables in the reduced array become Xn p., Xn T, Xn P 7f[ and p 0 . 

1 Xj 

This procedure is essentially the same as that described in Refs. 3 and 4. 

KINET 

This routine determines reaction rates 7 ascertains the controlling reac- 
tions, on the basis of which it forms the a ^ transformation matrix; prepared 
the reduced correction coefficients for the kinetic terms 7 and combines them 
with the previously determined coefficients. 

RE RAY 

This general purpose inversion routine inverts the reduced coefficient and 
multiplies the result times the set of errors yielding an unconstrained set of 
corrections for the reduced set of independent variables. 

CRECT 

The unconstrained corrections for the remaining independent variables are 
calculated by this routine. For species which are important in mass balances 
the logarithmic corrections are changed to linear corrections, i.e. 

APj = £j- A in Pj (113) 


29 


in order to emphasize the linearity of the mass balance equations; all cor- 
rections are checked with regard to the constraints of Eqs . (108) or (109) 

and the minimum multiplier is determined, and corrections are appropriately 
reduced and added to the previous trial balance. 

PROPS 

This routine evaluates transport properties according to equations of 
Ref. 1 as well as the derivatives appropriate to the boundary-layer solution 
(see Section 4) . 

These brief descriptions serve to give an overview of the routines in- 
volved in the state calculations. The program has been extremely successful 
on a wide range of problems and convergence has usually been quite satisfactory. 


SECTION 4 

THERMODYNAMIC PROPERTY EVALUATION 


Once solution is obtained for a given set of input parameters it is still 
necessary to evaluate a variety of properties which are state dependent. Some 
of these require simple summations and no particular discussion is required. 
Others involve either first or second derivatives of the state with respect 
to certain input variables. It is these properties with which the present 
section is primarily concerned. One of the incidental advantages of the 
Newton-Raphson procedure over other direct search or optimization (gradient) 
methods is that most of the information necessary for derivative determina- 
tion is already compiled when solution is achieved. 

Consider the closed system / assigned-enthalpy-pressure option. The re- 
duced set of independent variables used in the formulation are in p^, p^ # 
Xn(P^), and in T. In the program assignment, the independent variables are 
h, p, and a^, the last term being the gram-atoms of "element" i in a gram of 
system. The equations which implicitly describe this system were summarized 
in Table I. The problem reduces to a set of equations 

f i (x i ,Y i ) = 0 (114) 

where the might be considered as the reduced set of variables in p^, 

p lt in{PV\) , and in T, and the y^ are h, P and au. The goal is to achieve 
derivatives of the form (dx/dy) at constant all other y, for example 
(dT/dh) p , a ^. The general procedure for doing this is described in Appendix A. 
The appendix actually considers both first and second derivatives, although 
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in practice only the first derivatives have been calculated. For the first 
derivatives 


l 


i 



i 'H 


l* f k\ 

-l 

/ af k\ 

r x i . 



\ /V'Yi 


\ }x L . y A . 


i'^i 


i 'H 


(115) 


The matrix which must be inverted is identical with the error-coefficient 
array previously developed. The last vector term in the above equation is 
readily evaluated based on the error equations. Considering the error equa- 
tions of Table 1, the following table indicates the important derivatives for 
the assigned enthalpy, pressure and elemental composition problem. 


Derivative 
\ with re- 
\spect to 

Equations 

P 

h 

a k 

Enthalpy 

0 

-p m 

0 

Pressure 

-1 

0 

0 

Conservation 

0 

0 

-P57| 

Equilibrium 

0 

0 

0 


It should be noted that the reduction of independent variables that occurs 
when the equilibrium equations are substituted into the remaining equations 
does not affect these derivatives because P, h, and a^. do not appear 
in the equilibrium relations. Also the product is taken as a single 

independent variable (of the x-set) and thus does not enter into the partial 
derivative terms tabulated above. For the given example, derivatives of 
in T, 4n P 7/1, in p A , and p tf are determined with respect to h, P, and a T 


in p i , and p^ 

The derivatives of i n p. are obtained by applying 


k # 


c™> AH . 

d in Pj = 2_, v ji d in Pi + rt 3X11,11 (116) 

i 

Derivatives of composition-dependent properties can then be obtained by con- 
ventional procedures. It is perhaps worth noting the relation between some 
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of the more conventional thermodynamic properties and the base set of variables. 
Specifically 


3 l n T 

ah 


-l 

P#a x 


(117) 


a in gj 
a in P 


= p 


T/ ou 


fa in pjn\ 

a /n p yn\ 

a /n t\ 1 

a in T 

dP L 

(_l (h,a k 

l 8h UJ 


ap L „ 

' ,h ' a k-l 


- 1 


IL 

In 

M\ 

fa jen p^l 


a Jin T\ 


In 

T ) P,a k 

1 dh J 

!p. a k 

(ah ] 


-l 

P , cu 


(118) 

(119) 


Y = 


a in p 
9 in p 


a + 


a in 
d in 


m \ - -| a in a 1 

P >T,ou % L Xn i,a v J 


( 120 ) 


are useful relations. A convenient reference to other thermodynamic deriva- 
tives is listed under "Thermodynamic Formulas" in more recent editions of 
the Handbook of Chemistry and Physics. 

Similar relations to those given above can readily be generated for solu- 
tions at assigned entropy and pressure. For other cases it is often simplest 
to formulate the enthalpy, pressure, composition derivative array after obtain- 
ing solution based on a different set of constraints. This array is then used 
as indicated in the example above. 


SECTION 5 

SUMMARY AND CURRENT STATUS 

In the preceding sections a general chemical state procedure has been 
developed and mathematically applied to a number of open and closed thermo- 
dynamic systems. An effort has been made to provide a relatively general 
approach to the problems associated with such solutions and to indicate, via 
examples in some cases, means of circumventing them. A very brief discussion 
of the mechanics of the solution served to introduce the program and subprograms 
involved in the computer analysis. Some of these routines are quite general in 
their present formulation, others directed toward specific systems. 
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Currently all equilibrium aspects of the program are fully operational 
for general chemical systems. This includes the various closed and open sys- 
tem options, the shock wave relations, the surface coupled boundary layer mass 
balances, the bookkeeping involved with treating appearing and disappearing 
atomic elements, and the property and property derivative calculations. The 
KINET routine currently treats only the heterogeneous reactions associated 
with graphite oxidation and reduction. The generalization of this routine 
following the detailed approach presented in this report is a major recom- 
mendation of this report. 

The report has discussed in rather general fashion the treatment of gen- 
eral chemical systems. The ultimate program which should evolve from this 
study will be a General Nonequilibrium Ablation Thermochemistry (GNAT) program 
designed for treating the problems associated with equilibrium and nonequilib- 
rium at and above ablating surfaces. 
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TABLE I 


SUMMARY OF EQUATIONS AND RELEVANT INFLUENCE COEFFICIENTS* 

1. Gas phase Equilibrium: for all non-base gas phase species 

© 

ERROR = - in K pjt + Y j( ~ ^ 

i 

Derivatives with respect to; 


in p.: 5 j • , ( - 1.0 if j = j ' ; 


= 0 if 3 ¥ j') 


y . : — v • , ■ 

]'i 


or 


> © 


Pj : 0 


P i : ° 


<n P% 0 

-H it + Z v H 

in T: J — L-i-i 


RT 


2. Condensed Species Equilibrium: for all present non-base condensed 

© 


ERROR * - l 


n v, - E v i-i y i 


species 


Derivatives with respect to 


5n p. 


^i : 
{ or 
P i s 


P i: 


9. n V7f{\ 


9 n T: 


' i 

0 

0 

0 


>® 


— n * , +■ v m , . h , 
9 * 9 ' l j 


RT 


Notes in circles and indexing conventions are at the conclusion of the table., 


34 


TABLE I (continued) 


3. Surface Equilibrium: for the condensed species with the smallest 

algebraic error having a fail temperature 
equal to or greater than system temperature 

Z © 

v i*i^i if n on-base 

Xi ** 

1 

ERROR = - if i* is base species 

Derivatives with respect to: 


in py. 
/ 

Y i ! 
< or 

p i ! 

\ 

in Pft?: 

in T: 
in T: 


~ V i*i 


> © 


r 

h.4, + h 


\e* I i v -e*i H i 
RT 


if 

i* 

is 

non-base 

if 

i* 

is 

base species 


This constraint is deleted if the chosen species is a present condensed spe- 
cies, i.e., if the system temperature equals the species fail temperature. 


4. Closed System Enthalpy: 


ERROR 

= -Ptffh 

+ Z p j H i 




Derivatives with respect to 

lr\ p j : 

p j H i 


j « 

p.H. 
*1 l 



Hi j 


P* : 

H i 


in P$: 

-PJ?!h 


in T: 

Z c 

Pi P± 
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TABLE I (continued) 


5. Closed System Entropy: 
ERROR = -P^S + 




j»i,i 


Derivative with respect to 


In Pj! Pj (Sj - R) 


Y i : p i (S i “ R) 

or 

p i : S i 


p i : S i 

in (P#J : -P/7fs 

in T: p^C 


i» j 


6. Pressure: 


ERROR = -P + 


I 


3 ' 


Derivatives with respect to 


in PjS Pj 


y . : p . 

J X *1 


Pi* 0 


© 


y r - 
in (P #0 : 0 


in T: 0 
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TABLE I (continued) 


7 „ Oblique Shock Momentum: 


(P 


ERROR 


= - p i " 


U.COS0 ) 3 V -1 

\y ■ + L P j + <Pl u l co 80 l> 


a RT 

P57! 


Derivatives with respect to 


In Pj? Pj 


y±: Pi 

or 


© 


Pi? 0 


P i : 0 


o TR'P 

in P% - (p^cose^ 2 || 


jen T: (p.u COS0 ) 2 RT 




8. Oblique Shock Energy: 


r (u 1 c OS 0 1 ) s r 

ERROR = “P % |^h 1 + ^ J + ^ P jH j + -j (p 1 u 1 cos0 1 ) 




Derivatives with respect to 


X n p. : p.H. 

3 * 3 3 

( \ 

y p.H. 

2 1 *11 

or 




Pi : H. 


p je ! H i 


, ^ r (u.coso ) a 

In (P57!) s -p^l + J - j (p^^cos©^ 

£n T: (P^cose ^ 3 + T ]T c p .Pi 


a iRTl 2 

PWl 


i.j.i 1 


(RT ) 2 

pm 
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TABLE I (continued) 


9. Closed System Mass Balance: 
ERROR = -P 771CL. 


+ z 


V ji‘ P D 


Derivatives with respect to 


in p . : v . • , p . 


/ \ 

Y i ! 6 ii-Pi 

\ or / 

p i s 6 ii- 


p * : 

in (P»i) : -PJfau, 

in T: 0 


10. Surface Mass Balance Coupled to Boundary Layers 

ERROR = -P^ g a i , + (pv) (P^Y f + ^ p j v ji + p i' 


(pv) = m g + m 


E + E p i v <i' + V V E p j H j 

j.ig l 

v-E 


p iV p \ 


j£,i 


Derivatives with respect to 


Pj : 

(pv) 

Vji, + 

Si' + Y h i , H sJ 

Y i : 

or 

_(pv) 

6 ii- + 

x ii. * V*i] 

p i : 

1 - (\/P^g) 

F V» i . + Z 

j,1 c 

P^ ! 

v jei* 

- (^/P!7J ) 

p Vf,. + ! 


li’ 




+ (m + m ) y _ I P/71 + — 1 


In (P WJ: ^ _a i + (m g + m,) Y f J 

E p j H j 


g ~TW r 


PjVji, 




i n T: T 

n i 


3,1 
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TABLE I (continued) 


11 . 


Surface Mass Balance Coupled to Convective Coefficients: 


ERROR = -P 7 n g a z ^' + (pv) Pj v ji' 


l V li' + p i 1 


3*i r 


V p. IjjJl. 

^ 3 F^/P 


E p^j 

f = Ha 


E 

j,i g 


* . E 

m + m i , i 

(P V) = -3 c c 


P e U e C M E P-j^-i 

j .i. 


Derivatives with respect to 


in p j : 


(pv) + p* v . . , ,p. + 


Fj | ji'^j 


l V,.!. ^ 

L J TP* / 


L i'.t 


f*,/f 


(1 


F/F* ) 


V” i.i 

+ ^Pj.Vji.— 


E %p i i 




E p-ift 


3»i, 


3 J 


_ 2 £i 


£ 

j*i 


or 


(pv>+ p* 6 ii- p i' + 


E 


P D- 


L 3 M, 


F* ,/F 


(1 - 


F/P* ) 


L p j- v ji'T- 


3,1 


E P+*L 


3,lg 3 3j 3,lg 


E p* 


> © 


p. : 6 


1 11 


_ y D 

ii 1 P 7 h l — i P j v ji‘ 


3 ,i, 
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TABLE I (concluded) 


K 


p i : v l i' 


T p.7ft- r—} 
h* p ] ] i/i 


Z p j v ji’ 




in (P^ g ) : -^ g a z ., 

in T: 0 


INDICES 


i,i' 

j" 


i 


c 


E 




g 


base species or “elements" 
gas phase non-base species 
condensed non-base species 

included for condensed base species only 
included for gas phase base species only 
implies successive summation over all gas phase species 

NOTES 


1. The variable can have the following specifications 

= in p^, for gas phase base species 

= 0 , for present condensed base species or base species represent- 

ing a non-present element (in this case p^ becomes the 
system variable in place of y^) 

< 0 , for non-present condensed base species 

variable, for non-present base species representing a present element. 

In all but the second instance y i is a necessary but unknown variable 
of the system of equations. 

2. The variable p^ is used in lieu of y^ if the base species, i, is a 
present condensed species or represents a non-present element. 

3. Note the i if i* is a base species, likewise, = 6 -j_i 
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APPENDIX A 

DETERMINATION OF A SET OF FIRST AND SECOND DERIVATIVES 
IMPLICITLY DEFINED BY A SET OF NONLINEAR ALGEBRAIC EQUATIONS 


A set of K nonlinear algebraic equations in K variables and j 

parameters y ^ , can be expressed generally as 


f k (x i' Yj) = 0 


1 £ k £ K, l^i^K 


( 1 ) 


For any set of values for the parameters, y ^ , the ’ x^ are implicitly defined 
by the Eqs. (1). In many instances it is also desirable to obtained the ef- 
fect of changes in the parameters y_. on the variables x^, that is dx^/dy^ 
In the present case dx^/dy^dy^, are also sought. 

By differentiating Eq. (1) considering all x^ and y^ independent, 
the following partial derivatives can readily be obtained 


f k ' f k ' f k 


' f k . ' f k 


c. y. x.x.. y .y . x. y . 

i ii* *3 3 1 


The total differential of Eq. (1) is 


df = f dx + f dy = 0 

x i Yj 3 


( 2 ) 


where the repeated index implies summation. Dividing Eq. (2) by dy^ , , 
holding all other y^ constant, but allowing the x i to vary in accord with 
the Eqs . (1) , yields 


df, dx. 

h = f i + f 

dy . t x k dy . . r k 

* 3 ' x ± *3 Yj i 


= 0 


( 3 ) 


where a convention has been adopted whereby df,/dy. , implies (df v /dy. t ) 

K J x j y . 

j ^ j* but with x. varying to satisy Eqs. (1) and f ^ implies J 

l yj * 

(d f v /dy_. , ) „ v # j / j 1 . This set of K linear (in dx^dy^,) algebraic 


- . 3 

equatxons can be solved to yield 


dx . 

4 


( 4 ) 


X il 
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which represents the desired result for first derivatives . Taking the dif- 
ferential of Eq. (3) yields 



E iSL 

Vi- 9y :' 


dx . + f 

i* k 


x i y j 


dx. Sx. 

dy. + f. d - 5 —^- 

9y V 3 \ l Sy j' 


+ f,. 


y .y . 
D 


dy j 


+ f,. 


y . . x . , 

1 i' 


dx . . 

l 1 


= 0 


(5) 


If this relation is divided by dy.. at constant all other y. but with 

D D 

x. varying in accord with the Eqs . (1), there results 




9y j* ay j' 


= f, 


Sx. 3x. , 3x. d 3 x. 

i_ Li + f i_ + f _ .1 

: dy.. 3y .* k Sy.. k 3y.,3y .. 

X i X i' 3 3 X i y j* 3 x i 3 3* 


+ f,. 


+ £,. 


y j* Y D 


&X i' 

ay-s, 


= o 


i**-p yj- x i- 3i 


(6) 


Solving for d 3 x./dy.; , dy yields 
1 3 J* 


a 3 x.. 




X . J 


f 


\ x i x i- 3 


dx. 

r—~ + f v 

5y-i . X 


dx i( 


y r x i ./ ay 3* 


dx . 

+ f k ay - ? + f k 

x i y j* 3 *j' y j* 


(7) 


which represents the desired result for second derivatives. Note that the 
first term in the brace involves a triple sum (over k, i and i 1 ) and thus the 
time required to generate a single term is proportional to K a . To generate 
all terms in the three-dimensional array of second derivatives requires a time 
proportional to K 4 . 
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